Large-scale magnetic fields from hydromagnetic turbulence in the very early universe 
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We investigate hydromagnetic turbulence of primordial 
magnetic fields using magnetohydrodynamics (MHD) in an 
expanding universe. We present the basic, covariant MHD 
equations, find solutions for MHD waves in the early uni- 
verse, and investigate the equations numerically for random 
magnetic fields in two spatial dimensions. We find the forma- 
tion of magnetic structures at larger and larger scales as time 
goes on. In three dimensions we use a cascade (shell) model, 
that has been rather successful in the study of certain aspects 
of hydrodynamic turbulence. Using such a model we find that 
after 0(10^) times the initial time the scale of the magnetic 
field fluctuation (in the comoving frame) has increased by 4-5 
orders of magnitude as a consequence of an inverse cascade 
effect (i.e. transfer of energy from smaller to larger scales). 
Thus at large scales primordial magnetic fields are consider- 
ably stronger than expected from considerations which do not 
take into account the effects of MHD turbulence. 



I. INTRODUCTION 

It has been suggested that primordial magnetic fields 
might arise during the early cosmic phase transitions [Q , 
and recently it has been shown that magnetic fields are 
indeed a stable feature of the electroweak phase transi- 
tion . In a first order phase transition magnetic fields 
could also be generated when bubbles of the new vacuum 
collide, whence a ring of magnetic field may arise in the 
intersecting region j|] . It has also been suggested that at 
very high temperatures the ground state of a non-abelian 
particle theory is a "ferromagnetic" vacuum with a per- 
manent non-zero magnetic field The plasma of the 
early universe has a high conductivity so that a primor- 
dial magnetic field would be imprinted on the comoving 
plasma and would dissipate very slowly [H. Such a field 
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could then contribute to the seed field needed to under- 
stand the presently observed galactic magnetic fields [|j , 
which have been measured both in the Milky Way and 
in other spiral galaxies, including their halos. Typically 
the observed present day magnetic field is of the order of 
O{10-^) G. 

Locally the primordial field could be very large; it is 
limited only by the magnetic energy density contribution 
to primordial nucleosynthesis, which constraints the field 
at the time of nucleosynthesis by pB = B^^^^/8n < O.Spi,, 
where pi, is the energy density of 3 species of massless 
neutrinos M. Because flux conservation implies that 
Sims ^ R , where R is the scale factor, the field could 
have been much stronger at earlier times. On dimen- 
sional grounds, a typical value of the magnetic field fluc- 
tuation should be Brms ^ T^, so that at the time of 
the electroweak phase transition one could locally obtain 
fields as high as 10^** G. Depending on how such a strong, 
random magnetic field scales at large distances, it could 
1^ be the seed field needed to explain the magnetic fields 
observed on the scale of galaxies and larger. 

However, even assuming that a primordial magnetic 
field is created at some very early epoch, a number of is- 
sues remain to be worked out before one can say anything 
definite about the role of primordial fields in generating 
galactic magnetic fields. At earliest times magnetic fields 
are generated by particle physics processes with length 
scales typical to particle physics. If the inflation hypoth- 
esis proves correct, then after inflation rather long cor- 
relation lengths are possible The question is if it is 
at all possible for the small scale fluctuations to grow to 
large scales, and what exactly is the scaling behaviour 
of firms or the correlator {B{r + x)B{x)). Even in an 
inflationary scenario it would be of interest to see if the 
relatively large scale can grow even further. To study 
these problems one needs to consider the detailed evolu- 
tion of the magnetic fleld to account for such issues as 
what happens when uncorrelated field regions come into 
contact with each other during the course of the expan- 
sion of the universe. In general, turbulence is an essential 
feature of such phenomena. These questions can only be 
answered by considering magnetohydrodynamics (MHD) 
in an expanding universe. [|lO|| It is the main purpose of 
this paper to investigate the subsequent development of 
the primordial field. Expressed in a general way, our con- 
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elusion turns out to be that MHD turbulence is opera- 
tive, and hence the scale of magnetic fields is considerably 
larger than one would expect if MHD turbulence was ig- 
nored. This means that the previous estimates of the 
strength of the primordial magnetic field "today" need 
to be reconsidered. 

We begin by posing the basic equations and consider 
certain simplified models. A full H-3-dimensional nu- 
merical simulation would be desirable, but is beyond the 
scope of the present paper. In Sec. II we derive the rel- 
ativistic MHD equations for relativistic plasma, which is 
appropriate for the very early universe. In Sec. HI we 
discuss the appearance of waves in relativistic MHD. To 
illucidate the various MHD effects pertaining the early 
universe, we also present numerical solutions to the MHD 
equations for a two-dimensional slice. In Sec. IV we study 
a cascade model, that reflects important properties of 
fully three-dimensional turbulence. The cascade model 
has been rather successful in ordinary hydrodynamics. 
We find that in the early universe magnetic energy is 
transferred from small scales to large scales. We also 
compute the correlation function {B{r + x)B{x)) in the 
cascade model. In Sec. V we offer an interpretation of 
our results. 



II. RELATIVISTIC MHD IN THE EXPANDING 
UNIVERSE 

We begin by presenting a derivation of the fully gen- 
eral relativistic MHD equations (see also ref . |0 , where 
further references can be found), which we rewrite in a 
form suitable for our numerical work. We consider the 
early universe as consisting of ideal fluid with an equation 
of state of the form p = ^p, where p is pressure, p the 
energy density, and the speed of light is set to unity. We 
further assume that the fluid supports a (random) mag- 
netic field. The energy-momentum tensor is then given 

by 



T'^" = ip + p)^^ +pg' 



(1) 



and [/'^ is the four-velocity of the plasma, normalized as 



-1, and Fai, = dnA^, — d^A^ is the electro- 



magnetic field tensor. Note that, as long as diffusion can 
be neglected, the presence of the magnetic field does not 
change the equation of state. 

The magnetic energy is assumed to be much smaller 
than the radiation energy, so that it can be neglected as 
far as the expansion of the universe is concerned. We 
therefore assume a flat, isotropic and homogeneous uni- 
verse with a Robertson- Walker metric ds^ = —dt^ + 
R^{t)dx.^. Although the magnetic held generates local 
bulk motion, this may still be consistent with isotropy 
and homogeneity at sufficiently large scales, e.g. if the 
magnetic field is random, i.e. statistically homogeneous 



and isotropic on scales much larger than the intrinsic cor- 
relation scale of the field. Even very large magnetic fields, 
together with the ensuing very fast bulk motion, might 
not contradict isotropy and homogeneity. The equations 
of motion for the fluid arise from energy-momentum con- 
servation 
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The Maxwell equations read 
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We define in terms of the electric and magnetic fields 
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where latin letters go from 1 to 3. With this definition 
the expression for the total energy has no i?-factors and 
takes therefore the familiar form 
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where j — U'^. 

In order to solve (||) and (^) we rewrite the equations of 
motion explicitly in 3-1-1 dimensions. We start by writing 
(I) as 
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9"'-^-F"'9..r^ (6) 



where \/—g = R^, and the nonvanishing Christoffel sym- 

= iR/R)S) = r;„. It 

, because then the nor- 



bols are pO- = RRS,j and r„j 
is useful to define [/* ~ jR^^v^ 
malization U^U^ = — 1 gives the familiar form for the 



Lorentz factor 7 = (1 — v^) 
For p = i we obtain 



dt 



^1 



(V-v)(i?4S)- (v V)(i?4S) 
V(i?V) + (i?'J) X (i?2B)], 



(7) 



where S = (p -I- /3)7^v. It should be noticed that in this 
equation all quantities are scaled by the appropriate pow- 
ers of R. Thus, e.g. i?'*S is expected to be independent of 
i?, because p + p scales like 1 / R^, and v is expected to be 
independent of R. Also, V occurs always multiplied by 
or, alternatively, the operator d/dt is replaced by 
itself multiplied by i?, which means that time is replaced 
by conformal time t = J dt/ R. To emphasize this, it is 
convenient to introduce new scaled "tilde" -variables, 



S 
J 



and E = 



~P = R^P, 
i?2E. 



B = R^B, 



(8) 



It should be noticed that v is not scaled. Equation (|^) 
can then be written 
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8S . _ 

-^ = -(V-v)S-(vV)S- Vp- 

ot 



J X B. 



(9) 



For fi — we obtain, using scaled quantities 
47V dt dt 

J E 



r2H/ + ^+v-Vln(p7^) 



dt 
V-v = 



(10) 



In order to solve this equation numerically with an ex- 
plicit code we need to eliminate the time derivative 
d\nj^/dt. To this end we first solve the normalization 
condition for 7^, 



(11) 



where we have used p + p ^ |/9 for later convenience. We 




then differentiate 
dlnj^ ] 



dt 



dS^/di 7^ - 1 dh-ip 
■ di ■ 



72(272 -l)(p + p)2 ^2 



(12) 



Combining (0) and (|l|) we obtain a final equation suit- 
able for numerical work 



272 + 1 d In p 
472(272 - 1) di 



dS^/di 



(|P7) (27^ - 1) 
V • Vln(/572) - V • V 



J E 
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In this equation only scale invariant quantities enter. 
The Maxwell equations can be written explicitly as 



„ - 
— ^ = -V X E, 
dt 



V - B = 0, 



and 
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-V X B 



dH, 

'dt' 



V • E = Anpe 



(14) 



(15) 



where Pe is the charge density and jj^ 
E = -V X B, 



B?Pe- Further, 
(16) 



which is valid in the limit of high conductivity pi[ . 
Again, these equations have the natural scaling proper- 
ties with respect to powers of R. We emphasize that in 
the relativistic regime the displacement current, —dE/dt, 
cannot be neglected. However, in all cases considered we 

were able to solve for — E = v x B + v x B iteratively by 

evaluating v and B from the previous iteration. 

The equation of energy conservation is T^'^-i, = 0, or 

^§iiR'TOO) + ±.TO^+RRP^=0, (17) 



but since T"^ = 0, we have T^^ = T'j/R^ = -T°o/R^ 
T^^ /R^ and therefore the energy equation is 



or integrated over the whole space 

dR'^E tot 



dt 



0, 



where 



Etot = / T^'^d^x - /-^""^ 



(18) 



(19) 



(20) 



Hence R'^Etot is conserved. 

The conclusion from the above expressions is thus that 
the MHD equations in an expanding universe with zero 
curvature are the same as the relativistic MHD equa- 
tions in a non-expanding universe, provided the dynami- 
cal quantities are replaced by the scaled "tilde" variables, 
and provided conformal time t is used. The effect of this 
is, as usual, that the expansion slows down the rate of 
dynamical evolution. 

It should be noted that the velocity v is the bulk veloc- 
ity. Thus, in general we expect that v is nonrelativistic. 
This is physically reasonable since, although the gas par- 
ticles move with velocity near unity, we expect no strong 
collective effects which could give rise to a relativistic 
bulk velocity. The equations for nonrelativistic bulk mo- 
tions of a relativistic gas are given in Appendix A. 

In the early universe conductivity is large, and hence 
the diffusion length is also large. The conductivity of the 
isotropic relativistic electron gas, which interacts with 
heavy (non-relativistic) ions, is related to the Coulomb 
scattering cross section and reads |ll2| 



47r(Tcolin-e 



T 
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(21) 



where u>p is the plasma frequency, (TcoIi is the collision 
cross section and a is the fine structure constant. This 
result is valid for fields smaller than the critical field 
Be — ml/e = 4.41 x lO^'^ G, above which the elec- 
trons cannot be treated as free, and the conductivity ( |2l| ) 
should be multiplied by a factor B/Bc- On dimensional 
grounds, conductivity of the fully relativistic Standard 
Model gas will also scale as cr ~ T. The expansion rate 
of the radiation dominated universe is given by 



R 2t 



90 Mp 



(22) 



dx3 



where is the number of the effective degrees of free- 
dom, and Mp = 1.2 x 10^^ GeV is the Planck mass. 
Equation (E3) also provides the time-temperature rela- 
tionship |13||, and the inverse the length scale of the uni- 
verse. A measure of the importance of the diffusion is 
the magnetic Reynolds number, which may be defined as 
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Re = Lva, where L and v are respectively the typical 
length scale and velocity in the system under considera- 
tion. A Reynolds number less than one means that dif- 
fusion dominates. In the early universe, say at the elec- 
troweak phase transition Tew — O(IOO) GeV where in 
the Standard Model 5* ~ 106.75, the Reynolds number 
is huge; typically 



Reu ~ vaH 



Mp 



10' 



(23) 



where v has been arbitrarily chosen to be 10^^. In this 
sense the very early universe is almost a perfect conduc- 
tor. Also, the extremely large value of the Reynold's 
number indicates a turbulent situation, which we shall 
find by other methods later. 



We now seek a wave solution which, because of the 
structure of (|2^)-(|2^), must contain the scale factor R to 
the power —2, 



and 



(5B = — ^ exp i (k • X — uji) 
R-^ 



V = vo exp i(k • x — wi), 



const. ~ 
6p= ^4 expz(k ■ X - cjt). 



(27) 



(28) 



(29) 



where bp and Vg are constants. These expressions satisfy 
the basic fluctuation equations (p3)-(Efl) with 



III. ASPECTS OF RELATIVISTIC MHD 

A. Magnetohydrodynamics waves 

Let us begin by first presenting some general consid- 
erations. In the framework of relativistic MHD in an 
expanding universe, we can still discuss waves. Although 
the equations exhibited in the previous section are con- 
siderably more complicated than their non-relativistic 
counterparts, the MHD-waves are linear perturbations of 
the standard cosmological background. Thus, the bulk 
velocity v must necessarily be small relative to the ve- 
locity of light. It therefore follows that the displacement 
current can be ignored. The background is homogeneous, 
and we assume the relativistic relation p = p/i for the 
background as well as for the fluctuations. The continu- 
ity equation, i.e. (p^), gives to the lowest order the well 
known result p = const /i?^. To the next order we get 



t 



dRHp_ 

at 



= 0, 



(24) 



where Sp is the fluctuation in p. Also, from (^ we get to 
lowest order in the fluctuations 

= -| ^R'ilsp + BSB) + l(i?2B)V(i?^5B). 

(25) 

Here 6S = |/0v, and (5B is the fluctuation of the back- 
ground field B, which is assumed to behave like ~ 1/ R^ . 
Of course, ^B is expected to have a similar scaling be- 
havior as a function of time, but it also has a spatial 
dependence. Finally, we have the fluctuation equation 



dR^SB 
dt 



V X (v X i?2B), 



(26) 



which follows from (|l4|), since the displacement current 
can be ignored for small bulk velocities. 



dt' 



(30) 




where <o is the initial time, in accordance with the results 
obtained in the previous section. 

We therefore see that with the scaling properties 
mentioned above the equations are similar to the non- 
relativistic case (provided we use the time i in (|3^) . Thus 
we find the group velocity du/dk — B/ -|- p. Because 
the scaling properties of B and ^/p + p with respect to the 
expansion of the universe are the same, it follows that the 
group velocity is independent of R. As far as the phase 
velocities are concerned, the same is true. Assuming the 
background field to be in the x-direction, then 5B and 
V are in the z-direction, as in the case of non-relativistic 
waves One then finds that the velocities are given 

by 



(31) 



Of course, these velocities are given in terms of the con- 
formal time t. If we express the result in the standard 
time t, then (^Tj) should be multiplied by \/R. It should 
be noted that the assumption of small bulk velocities can 
only be maintained if |B| ^ ^J]). If this condition is not 
satisfied, we cannot expect the nonlinear effects to be 
small. 



B. Two-dimensional slice 

Ideally, we would like to solve the MHD equations in 
three (plus one) dimensions. However, as indicated in 
Sect. II, this is a major computational task. We restrict 
ourselves therefore to a two-dimensional slice only. The 
main conclusion will be that much of the qualitative be- 
havior of nonrelativistic MHD carries over to the case of 
relativistic MHD. 

We solve (|) and (|l|)~(|l6|) numerically using 6th order 
centered differences to compute the spatial derivatives 
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and a 3rd order Runge-Kutta scheme for the timestep. 
We adopt random initial conditions for B. In order to 
guarantee that V • B = for all times, we advance the z- 
component of the vector potential by means of the equa- 
tion dAz/dt = • (v X B), where B = V x (A^ez). 
The initial is computed by solving V^A^ = —AnJz- 
Initially {t — to, i.e. t = Q) we put p to unity. Pe- 
riodic boundary conditions are adopted in the x- and 
y-directions. 

Our new equations (^) -f ([l^)--(|l^) are scale-invariant, 
so it is sufficient to solve them on a computational do- 
main with size L — \. The results for a different do- 
main size L' are the same, but taken at a different time 
t' = {L'/L)l 

Like in all turbulence calculations there has to be some 
diffusion to prevent the accumulation of energy at the 
smallest scale. In order to restrict the effects of dissi- 
pation only to the largest possible wave numbers we use 
hyperdiffusion, i.e. instead of the usual diffusion operator 
V^, we use an operator of the form —V* for the evolution 
of all variables. This technique is well known in turbu- 
lence research (see e.g. ||l^). Also, since this procedure 
is merely of computational relevance, we did not use the 
relativistic expressions. 

The minimum diffusion coefficient v we can afford 
is given by the empirical constraint that the "mesh 
Reynolds number" i?emosh — USx/v must not exceed the 
value 5 — 10. Here, Sx is the mesh size and J7 is a typi- 
cal velocity that includes the velocity of waves and bulk 
motions. As was pointed out before, in the early uni- 
verse the Reynolds number is very large, which means 
that the magnetic diffusivity rj = An/a should be much 
smaller than the adopted value of ly. In other words, in 
order to have realistic values of 6x has to be extremely 
small. However, the maximum number of meshpoints, 
N = L/6x, is limited by computer memory and time. 
Our present, rather exploratory calculations were carried 
out on a workstation, and so we restricted ourselves to 
A'max = 128. Even on larger computers we would never 
reach realistic values. This demonstrates the difficulty 
of a realistic simulation. It is obvious that numerical 
simulations with a low Reynolds number cannot provide 
a realistic picture of the early universe MHD. However, 
we believe they are useful in illustrating the qualitative 
features of the problem. 

The evolution of the magnetic field is compared in 
Fig. for lower and higher resolution. As time goes on, 
the coalescence of magnetic structures leads to the grad- 
ual formation of larger and larger scales. In the higher 
resolution case there are more small scale structures, but 
also here the development of large scale fields is evident. 
In turbulence research this phenomenon is known as an 
inverse cascade. Such a cascade processes are linked to 
certain conservation properties that the basic equations 
obey. For further details see ref. |T^. We mention here 
only a few important aspects. An inverse cascade ex- 
ists both in two-dimensional and in three-dimensional 
MHD turbulence. The only difference is that in the two- 



dimensional case it is an inverse cascade of the magnetic 
potential, whereas in the three-dimensional case it is an 
inverse cascade of the magnetic helicity density A • B. 
In fact, the conserved quantities in the two cases are 
j cPxA^ and / c?'^xA-B, respectively. For comparison we 
also mention that the difference between two and three- 
dimensional hydrodynamic (non-magnetic) turbulence is 
more drastic. In two-dimensional hydrodynamics there 
is an inverse energy cascade associated with the conser- 
vation of enstrophy (mean squared vorticity), which has 
no counterpart in three-dimensional hydrodynamics. 
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FIG. 1. Left column; magnetic field lines at different 
times at low resolution (64 x 64 meshpoints). Right column: 
magnetic field lines at different times at higher resolution 
(128 X 128 meshpoints). 

The significance of an inverse cascade is that it leads to 
a transfer of magnetic energy to larger and larger scales. 
This process is due to the nonlinear terms giving rise to 
a mode interactions. Energy spreads over different scales 
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until some balance is achieved where the kinetic and mag- 
netic energy spectra have a certain slope. In the ordinary 
MHD turbulence a relevant energy spectrum could be the 
Iroshnikov-Kraichnan spectrum |17|| , where the spectral 
energy goes like fc~^/^, or a Kolmogorov type spectrum 
like k~^/^. These different spectra describe equilibrium 
situations, but in any case it is clear that the spectrum 
will be very different from white noise, which has a fc+'^ 
power spectrum (see Sec. IV. B below) . The possibility of 
energy transfer from small to large scales via an inverse 
cascade could be of major importance in cosmology. It 
could provide a seed field at the parsec or kiloparsec scale, 
albeit at small amplitude. 

IV. A CASCADE MODEL 
A. Description 

The ultimate goal is to solve the basic MHD equations 
in three dimensions at high resolution, using random ini- 
tial conditions. Although we would be unable to cover a 
realistically large range of length scales, it is important 
to know whether dynamo action could be possible in a 
relativistic flow. This is a major task, which would go 
beyond the scope of this paper. To see the difficulties in- 
volved in such a program, the reader should recall the dif- 
ficulties in making long term weather predictions based 
on the Navier-Stokes equations. Therefore, in order to 
demonstrate some of the anticipated behavior of the full 
1+3 MHD universe, we now study a cascade model of 
hydromagnetic turbulence. 

In ordinary hydrodynamics and hydromagnetics many 
properties of turbulence, in particular those related to 
energy transfer and to the spectral properties, including 
small intermittency corrections, have been studied suc- 
cessfully using a simple cascade model Q . This is true 
not only qualitatively, but also quantitatively, which is 
the reason why the cascade model is now much used in 
studies of nonlinear physics (see e.g. and references 
therein) . 

The basic idea is that the interactions due to the non- 
linear terms in the MHD equations are local in wavcnum- 
ber space. In fc-space the quadratic nonlinear terms, such 
as V X (v X B), V • Vv, and J x B, become a convolution 
and have the general form ||2^ 

Nk(v,B) = J Ca,pVa,ip)Bp{p-k)d^p . (32) 

(There are similar terms also for the other two nonlinear- 
ities.) where Caf3 = CQ^(k) is a tensor which is linear in 
k. Interactions in fc-space involving triangles with similar 
side lengths have the largest contribution, as discussed 
in [Q. This has led to the shell model (see e.g. p9| ] 
and references therein) , which is formulated in the space 
of the modulus of the wave numbers. This space is ap- 
proximated by N shells, where each shell consists of wave 



numbers with 2" < fc < 2"+-'^ (in the appropriate units). 
The Fourier transform of the velocity over a length scale 
fc~^ (fc„ = 2") is given by the complex quantity Vn, and 
Bn denotes a similar quantity for the _B-field. Further- 
more, the convolution is approximated by a sum over the 
nearest and the next nearest neighbors, 

2 

Nn{v,B)^ ^ CtjVn+iBn+j. (33) 
2J = -2 

Here v and B have lost their vectorial character, which 
reflects the fact that this model is not supposed to be an 
approximation of the original equations, but it should be 
considered as a toy model that has similar conservation 
properties as the original equations. Thus e.g. the en- 
ergy flow should be represented by these equations. It 
is quite remarkable that such models show several real- 
istic features, including intermittency corrections to the 
structure function exponents, and are therefore rather 
popular both in the absence ||l^ and in the presence pl| ] 
of magnetic fields. Therefore we propose to apply such a 
model also to the early universe. 

Velocity and magnetic fields are thus represented by 
a scalar at the discrete wave numbers fc„ = 2" (n = 
l,...,iV), i.e. fc„ increases exponentially. Therefore such 
a model can cover a large range of length scales (typically 
up to ten orders of magnitude) . The important conserved 
quantity is EtdR'^, where Etot = J T'^^d^x is the total 
energy. Using that the bulk velocity is nonrelativistic, we 
have 7 ^ 1 , so we can expand 7^ « 1 -f . Hence 




Since we are here mostly interested in the evolution of 
the magnetic field we ignore the detailed evolution of p 
and assume p « pqR~^. Thus, we require that 

J {^pov'^ + ^B^R^) d^x ^ const. (35) 

We use b = Bi?^ and construct equations for Vn and 6„ 
such that 

n— 1 n— 1 

In computing the conservation of the energy, the complex 
conjugate of this equation should be added. However, it 
turns out that the "complex energy" (exhibited in the 
above equation) is conserved by the following construc- 
tion. 

As pointed out, the main idea of the cascade model 
is to construct a set of equations that share the same 
basic conservation properties of the nonlinear (quadratic) 
terms as the original equations. Thus we write equations 
which mimic equations (^ and (p^), 

|po% = iV„(z;,6), (37) 
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— Mn{v b) (38) energy in fc-space. We used N = 30, which covers 

dt a range of length scales of approximately ten orders of 

^j^^^^ magnitude. As one can clearly see, magnetic energy is 

transferred from small scales to large scales. This is called 

2Nn{v,b) = jA:„(A + C)(w*+iU*+2 - inverse cascade effect. Such an effect is found in 

+ik„{B - |C)«„i<+i - (39) many nonlinear systems, for example in two dimensional 

-ikni^B + |a)«_2<-i - ^*-2^*-i): turbulcncc, relevant e.g. for the atmosphere. 



Mn{v,b) = ikn{A - C)«+i6;+2 - 6;+i<+2) 

+*/c„(B + iC)«_i6*+i - 6*_i<+i) (40) 

-ifc„(ii? - |^)«_2fe;-i - K-2<-i), 

with A, B, and C being free parameters. It is straightfor- 
ward to verify that "iJ^Vn^n + Y^ b'^Mn = 0, using that 
kn = 2". The i differentiations in Eqs. (|3^) and (^ are 
included to mimic closely the nonrelativistic form of eqs. 
(9) and (14). In the actual computations we have re- 
stored magnetic and kinematic diffusion terms, — i^fc^u„ 
and —rjk^bn, on the right hand sides of ( |37| ) and (^), re- 
spectively. We chose v = rj and as time goes on, lowered 
V gradually using the formula v > {Y, fc^l^nH^/V^tLx, 
where /cmax = 2^. This formula estimates the mini- 
mum amount of diffusion necessary to prevent the built- 
up of energy at the smallest resolved scale. We use a 

3rd order timestep, which is calculated via the formula 
,5t = 0.25min[(^fc2|y^|2)-i/2]_ 

B. Results 

The numerical study of the cascade model requires of 
course that the parameters A, B, C are fixed. This prob- 
lem turns out to be quite interesting, since it allows one to 
associate the cascade model with a dimension. In hydro- 
dynamics the parameters are fixed by taking into account 
a conservation law which is non-trivial in the dimension 
considered. In two dimensions, for example, one uses the 
requirement that the enstrophy is conserved, whereas in 
three dimensions the helicity should be conserved. In 
three dimensions Jensen et al. [|l9| used the values ^ = 1, 
B = —1/2, and C = 0, which we have adopted also in 
several models presented here. We compare the results 
with another set of parameter for which the quantity 

is conserved, in addition to the total energy |2^. This 
requires that A = 7/5, B = -1/10, and C = 1. The 
quantity Hm resembles the magnetic helicity A-B, which 
is important, because associated with it is the inverse 
cascade of magnetic helicity and energy |22j. In Fig. || 
we plot the spectral magnetic energy density EM{kn) = 
|6„p/(fci?^) computed for these values, with the initial 
field taken to be random (i.e., EM(k) = fc^). The reason 
we interpret this expression as the magnetic energy is 
that we know that ^ b'^bn enters in the conserved energy. 
However, Y ^ / = S dk/{k log 2), so EM{kn) is 
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FIG. 2. Spectra of the magnetic energy at different times. 
The straight dotted-dashed line gives the initial condition 
{to = 1), the solid line gives the final time (t = 3 x 10 ), 
and the dotted curves are for intermediate times (in uniform 
intervals of Alog(f - to) = 0.6). A = 1, B ^ -1/2, and 
C = 0. 

The quantity of paramount interest is the the magnetic 
field correlation function 

CB{r) = {B{r + x)B{x)), (42) 

which is related to the power spectrum via a Fourier 
transform, Csir) = J EM{k) cos{kr)dk. It is difficult in 
general to compute this quantity, due to the fluctuations 
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in the spectrum EMik). Therefore we have computed 
Csir) from the spectra of the cascade model by inter- 
polating EM{k) on a uniformly spaced mesh. This of 
course introduces some uncertainty. The result is shown 
in Fig. ^. Note the clear increase of the widths of the 
correlation functions. 




-0.4 I , ! 0.0 I , 

2x10^8 4x10^8 l.Qxlo"'^ 2.0x10 




,-5 



l.QxlO * 2.0x10 ■* 5x10 IxlQ * 

r r 

FIG. 3. Left column: The correlation function of B for 
three different times. Right column: The RMS magnetic field 
as a function of distance for three different times. 

Note also the anticorrelation at larger length scales. 
For a magnetic field this is natural, because if one con- 
siders the field in some region from a point far away from 
this region, the magnetic field in the region appears to 
be approximately a dipole. A negative correlation then 
arises because the field loop has to close. This would 
then basically be a consequence of div B=0, which has a 
tendency to lead to negative correlations. However, the 
cascade model of course has the difficulty that it does not 
really operate in ordinary space, but instead it is formu- 
lated in the modulus of k-space. Hence we cannot really 
investigate to what extent the condition div B=0 is sat- 
isfied, in contrast to the two dimensional case discussed 
in Sec. III.B. 

Another quantity of interest is the average magnetic 
field as a function of distance, 

B(r) = -1 y d^a:B(x), (43) 

where the integration is over a volume of size in D- 
dimensions. From this definition we have 



(BM') = ^ jd^'xj d^2;(B(x)B(y)), (44) 

where both integrations are over a volume of size . 
Thus the root mean square magnetic field 

i?. = (B(r)2)i (45) 

can be computed directly from the correlation function 
CB{r) via 

For a random field, Br behaves like r^"^, so the inter- 
esting question is whether this initial behavior changes 
as time passes. In Fig. |^ we show the results. There 
is a clear broadening of _B,. towards larger distances as 
time passes, as we would expect from the inverse cascade 
behavior. 

The determination of the width of the correlation func- 
tion above is not very accurate because of the fit involved 
in computing the Fourier transform of the spectrum. We 
shall therefore now introduce another length scale that is 
easier to compute, but whose value is similar to the width 
of the correlation function. The relevant length scale in 
turbulence theory is the so-called integral scale, which is 
the characteristic length associated with the large ener- 
getic eddies of turbulence. Roughly speaking one could 
view it as a measure of the coherence of the magnetic 
field, too. It is defined by 

^0== J 2Trk-^EM(k)dk J EM{k)dk, (47) 

which, in our cascade model, corresponds to — 
X]27rfc~^|6„P/^ |6„p. If the spectrum is random we 
get ^0 — |27rfc~;^^, where 27r/fcniax is the shortest length 
scale present in the model. This length scale in the ini- 
tial random spectrum is determined by the mechanism 
generating the primordial field. In Fig. ^ we show the 
evolution of Zo(t) in two cases, namely for the hydromag- 
netic A, B, C (circular points) and for the MHD A, B, C 
(diamond-shaped points). Although the two sets of val- 
ues for A, B, C do not yield identical results, we see that 
the curves are qualitatively similar. In both cases in- 
creases rapidly by 4-5 orders of magnitude, and there is 
a plateau structure. The MHD result (diamond-shaped 
points) has a plateau stretching to t/to = 10®, but for 
larger times keeps increasing. The increase of Iq by 
almost 5 orders of magnitude is important, because it 
could lead to magnetic fields at the present time at length 
scales comparable to one parsec. If we take the elec- 
troweak phase transition as the initial state, the QCD 
phase transition occurs approximately for t/to=10^. The 
maximum time t/to= 10^ reached in our simulation cor- 
responds to a temperature of 3 MeV, which is close to 
the nucleosynthesis. It should be noticed that from these 
results one cannot, of course, say anything about what 
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happens at later times. Therefore it could be that lQ{t) 
increases further, either by reaching new plateau(s), or 
otherwise. 
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FIG. 4. Evolution of the integral scale. The plot symbols 
denote uniform intervals of Alog(t — to) = 0.6. The circu- 
lar points corresponds to the hydrodynamic values A = 1, 
B = —1/2, and C = 0, the diamonds correspond to the MIfD 
values A = 7/5, B = —1/10, and C — 1, and the triangles are 
the widths computed for the correlation function. 

We also measured the integral scale of the magnetic 
field in the two-dimensional model of Sec. III.B and found 
a clear increase with time. For the three times plotted in 
Fig. 1 we found for the 64 x 64 case the values = 0.09, 
0.50, and 0.95, whereas in the 128 x 128 case values were 
^0 = 0.04, 0.28, and 0.76. The initial difference of a 
factor two is due to the different resolution. At later 
times the integral scales for low and high resolution are 
more similar. 

The evolution of Iq is not straight, but if we make a 
linear fit through the values given by the diamonds in 
Fig. U we find (ignoring the steep initial increase) 

lo{t)^ro{t/toy/\ (48) 

For further applications of the rough fit formula ( ^8|) it 
may be more sensible to express time in terms of tem- 
perature, T (X t~^/'^, so 

Zo(T)«ro(T/ro)-i/2, (49) 

where rg « 10~^; see Fig. ^. If we were to extrapolate 
to the present time we first have to fix the scale rg by 
physical arguments. The various models presented in the 
literature (lj|^ give characteristic scales for the primordial 
field when it is generated. This scale should be identified 
with the lowest scale in our calculations which, in the 
case of the shell model, is about 10~®. The scale tq is 
typically somewhat larger (10~^ in the shell model). The 
reason for this is presumably that a purely random initial 
condition is not consistent with the MHD equations. 

In order to clarify these point we take an example. If 
we assume that at the time of the electroweak phase tran- 
sition {T = 100 GeV) To was 10"'^ cm (the horizon scale 
was ~ 4 cm) then, using our extrapolation ( ^9[ ) , we arrive 



at a scale of 2 parsec. If we assume that the initial mag- 
netic field was 10^^ G, then the present day value would 
be 10^^^ G. Such values would lead to sufficiently strong 
seed magnetic fields to explain the field even in high red- 
shift galaxies by dynamo action Q. This extrapolation 
may be too naive, because the nature of turbulence will 
change as the universe cools down. Furthermore, at later 
times, when structure formation begins, gravitational en- 
ergy may lead to additional stirring and enhancement of 
turbulence in localized regions. 

V. DISCUSSION 

In the two-dimensional case we found that, starting 
from a small scale magnetic field, magnetic structures de- 
velop at progressively larger scales. This process of self- 
organization corresponds to an inverse cascade of mag- 
netic energy and helicity. Using then a cascade (or shell) 
model to study three-dimensional MHD turbulence we 
were able to follow this inverse cascade over much longer 
times. Such a cascade model has been rather successful 
in the study of hydrodynamic turbulence. 

The possibility of an inverse cascade means that the 
scale of fiuctuations of the primordial magnetic field in- 
creases much beyond its original scale given by particle 
physics. Taking the parameter as a measure of the 
coherence length of the magnetic field, we see that there 
is an increase in the coherence of 4 to 5 orders of mag- 
nitude. This means that previous estimates of the field 
strengths in various mechanisms for generating a primor- 
dial field should be revised accordingly. For example, let 
us consider the estimate by Vachaspati in ref . j^] . Taking 
the area average one has the estimate 

Br ^ gT^m, (50) 

where N is the number of steps needed to reach a given 
scale in terms of the "fundamental" scale at which the 
field is generated. In this case the fundamental scale is 
the electroweak scale ||]. Proceeding as in ref. [|l[ one 
has N ^ 10^^ today, if the relevant scale is of order 100 
kpc. However, due to the MHD corrections N is, from a 
conservative point of view, reduced because of the effect 
of turbulence, and one would instead have N < 10^^, 
which reduces the stochastic decrease of Br- It should 
be emphasized that from our calculations one can only 
say what happens up to a time of order 10^ ^eWi so 
presumably N is considerably below 10^^ today. 

The turbulent nature of the magnetic field may have 
interesting effects on the various phase transitions in the 
early universe. Also, the inherent shift of energy from 
small to large scales may be of interest in connection with 
the density fiuctuations due to the magnetic energy. 

Of course, the cascade model is a model of the real 
1+3 dimensional MHD turbulence. Its successful appli- 
cation in many, widely different nonlinear physical prob- 
lems suggests, however, that it might also be applicable 
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to the primordial magnetic fields of the early universe. 
Therefore we believe that its indication of the strong in- 
crease in the coherence scale of the primordial field should 
be taken seriously, and that further, more detailed stud- 
ies are warranted. 
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APPENDIX A: 

We give here the equations for the case where the bulk 
velocity is small. (The gas remains still relativistic, which 
is important for the scaling properties of p and p.) 



dlnp 
di 



-(v • Vlnp + V • v) 



J E 

p 



(Al) 



- = -v( — ^ + V-v) --Vln/?+^^, (A2) 



Dt 



Dt 



where D/Dt = d/di + v ■ V is the total derivative, and 

-^-Vx(vxB), VxB. (A3) 

dt 47r 
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